MODEL  BASED  CLASSIFICATION  USING  MULTI-PING  DATA 


Christopher  P.  Carbone  Steven  M.  Kay 

Naval  Undersea  Warfare  Center  University  of  Rhode  Island 

Division  Newport  Department  of  Electrical  and 

1176  Howell  St.  Computer  Engineering 

Newport,  RI  02841  4  East  Alumni  Ave.,  Kingston,  RI  02881 

(email:  carbonecp@npt.nuwc.navy.mil)  (email:  kay@ele.uri.edu) 


ABSTRACT 

This  paper  proposes  a  method  of  target  classification  using 
three  dimensional  (3-D)  data.  The  data  consists  of  multi¬ 
ple  realizations  (pings)  of  range  versus  bearing  plots,  so  the 
three  dimensions  of  the  data  are  range,  bearing  and  time  (or 
pings).  The  data  is  assumed  to  consist  of  independent  non- 
identically  distributed  complex  gaussian  noise,  and  a  target. 
The  Target  (TGT)  is  of  known  constant  size  (extent  in  range 
and  bearing)  and  known  speed.  The  TGT  power,  and  head¬ 
ing  are  unknown.  In  the  derivation  of  the  classifier  a  nor¬ 
malization  step  is  necessary  and  we  propose  an  approach  to 
the  normalization  of  multidimensional  (m-D)  data.  This  pa¬ 
per  contains  the  derivation  of  the  classifier,  a  description  of 
the  normalizer,  a  description  of  the  algorithm  that  follows 
from  the  classifier  and  simulation  results. 

1.  INTRODUCTION 

In  the  past  SONAR  data  has  been  processed  by  breaking 
up  the  3-D  data  into  one  (and  sometimes  2-D)  pieces  (see 
[1]  and  [2]).  This  has  been  done  to  simplify  processing 
and  reduce  computational  load.  The  normalization  step  has 
been  done  by  windowing  the  data  in  one  or  sometimes  two 
dimensions  as  well.  This  has  the  disadvantage  of  not  us¬ 
ing  all  available  information  to  perform  the  normalization. 
This  paper  will  present  an  algorithm  developed  by  using  the 
generalized  likelihood  ratio  test  (GLRT)  (see  [3]).  The  re¬ 
sulting  algorithm  operates  on  3-D  data  without  windowing 
or  breaking  the  data  up  into  one  or  two  dimensional  parts. 
This  algorithm  can  be  used  to  normalize,  cluster  (group  like 
threshold  crossings)  and  classify  the  data.  In  particular,  this 
algorithm  inputs  pings  of  data  that  have  been  beamformed 
and  outputs  a  list  of  possible  TGTs. 

2.  PROBLEM  STATEMENT 

Sonar  data  can  be  thought  of  as  existing  in  three  dimen¬ 
sions.  In  this  case,  the  3-D  data  is  formed  by  stacking  N 


range-bearing  plots  (see  Figure  1  with  N  =  3).  Each  range¬ 
bearing  plot  is  referred  to  as  a  ping’s  worth  of  data.  Other 
dimensions  could  be  added  such  as  doppler  but  in  this  pa¬ 
per  we  will  be  considering  the  discrete  3-D  case  of  range, 
bearing  and  ping. 


Ping  1 


5  10  15  20  25  30 


Beams 

Fig.  1.  Three  range  bearing  plots 

The  following  derivation  and  simulation  are  indepen¬ 
dent  of  units.  Throughout  the  paper,  range  is  in  lines,  bear¬ 
ing  is  in  beams,  and  time  is  in  pings.  Lines  could  be  any 
distance.  Beams  could  be  any  number  of  degrees,  and  the 
time  between  pings  is  not  specified. 
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3.  DERIVATION  OF  THE  ALGORITHM 


•  bsi  is  the  bearing  speed 


3.1.  Derivation  of  the  Classifier 

Consider  multiple  range-bearing  plots  (multiple  pings  of  range¬ 
bearing  data).  We  are  assuming  the  noise  is  independent 
non-identically  distributed  complex  gaussian  noise  (CN).  A 
range-bearing-ping  cell  of  CN  is  denoted  by  re  (mo,  no,  to)  = 
u(mo,  no,  to)  +jv(mo,  no,  to)  where  u  and  v  are  real,  inde¬ 
pendent  and  distributed  as  u  ~  N(0,  cr^/2),  v  ~  N(0,  cr^/ 2) 
and  a%j  is  the  power  of  the  cell. 

In  the  absence  of  a  TGT  we  assume,  the  “signal”  is  noise 
only, 

x(m,  n,  t)  =  w(m,  n,  t)  (1) 

where  w(m,  n,  t)  ~  CN( 0,  Pm,n,t)  f°r  ni-lines,  mm  0, . . . ,  M- 
1,  for  n-beams,  n  =  0, . . . ,  N  —  1,  and  for,  t-pings,  t  = 
0,...,T-1. 

When  a  TGT  is  present  we  observe 

x(m,  n,  t)  =  Si(m,  n,  t)  +  w(m,  n,  t)  (2) 

where  i  =  1,  2, ...  k  with  k  being  the  number  of  TGT  Mod¬ 
els  and  with, 


CN(0,  Ptgt)  (m,  n,  t)  e  A(m0,n0)i 

otherwise 

(3) 

Ptgt  is  the  TGT  power,  an  unknown  constant  and 

A(m0,  n0)i  =  {a(t  =  0);,  a(t  =  1  );,  ...a(t  =  T-  1)0 


The  TGT  speed  is  assumed  known,  but  the  heading  is  not, 
so  that  rsi  and  bst  are  constrained  but  not  known. 


i(t  =  O)*  = 


a(t  =  1  )i  = 


[mo  +  Or  Si,  mo  +  re  +  Ors^]  x 
[no  +  0  bsi,no  +  be  +  0  bst] 

[mo  +  lrsi,  mo  +  re  +  lrs^]  x 
[no  +  Ibsi,  no  +  be  +  lbs  A 


To  set  up  the  classifier  we  have: 

Ho  :  x(m,n,t)  =  w(m,n,t) 

Hi  :  x(m,n,t)  =  Si(m,  n,  t)  +  w(m,  n,  t)  (5) 

Assumptions  in  the  classifier: 

•  All  range-bearing-ping  cells  are  independent  (but  not 
identically  distributed)  random  variables. 

•  Pmnt  is  known. 

•  Ptgt  is  known  and  constant 

•  One  of  the  TGT  models  matches  the  true  TGT  head¬ 
ing.  See  section  5  on  relaxing  this  requirement. 

Under  Hi  :  p(x;  m0,  n0,  ,  Hi)  = 


nnn  0 

{m,n,t)^A{mo,no)i 

with  Smrnt  —  Ptgt  T~  Pmnt  • 
Under  Ho  :  p(x;  Ho)  = 

nnn 

(  4.\r-  a  t  \  '  r  mnt 

(m,n,t)eA{mo,no)i 

nnn 

(m,n,t)£A(mo,no)i 


—  | x(m,  n,  t)  | 

Pmnt 

—  | x(m,  n,  t)  | 

P mnt 


—  | x(m,  n,  t)  | 

P mnt 

—  \x(m,  n,  t)  | 

Pmnt 


a(t  =  T-l)i  =  [m0  +  {T  -  l)rsi, 

mo  +  re  +  (T  —  l)rs^]  x 
[n0  +  (T  —  l)6si, 
no  +  fre  +  (T  —  l)bsi] 


where, 


•  mo  is  the  start  of  the  non-zero  range  cells 

•  no  is  the  start  of  the  non-zero  bearing  cells 

•  re  is  the  range  extent 

•  be  is  the  bearing  extent 

•  rsi  is  the  range  speed 


As  a  result  we  have, 


p(x]m0,no,Si,Hi )  _ 


nnn  1 

1 1  1 1  1  l(m,n,t)eA(m0,n0)i  TrSmnt  C  ^  \  5mnt 

nnn  i  f 

II  II  ll(m,n,t)eA(m0,rao)i  tt Pmnt  ^  V  Pmnt 


Simplifying  we  have  = 


-4 exp  (EEE 


(m,n,t)£A(mo,no) 


Sex p  EEE 


(m,n,t)GA(mo,no 


Pmnt 

\x(m,n,t)\: 

P mnt 


with  ^ = n  n  n 

and  S  =  n  I!  Ilf, 


l(m,n,t)GA(mo,no)i  T^Smnt 

1 

(m,n,t)£A(mo,no)i  irPmnt 


Taking  In  of  both  sides  we  have  In  )  = 


ZIZm  \x(m,n,t) \2(-^ - tA  )  (10) 

/  ^x  ^  ,  /  N  \z  rant  ^mnt  . 

(m,n,t)£A(mo,no)i 


3.2.  Description  of  the  Normalizer 

The  range-bearing-time  data  is  normalized  by  estimating 
the  background  power  with  a  3-D  minimum  variance  spec¬ 
tral  estimator  (MVSE)  and  then  dividing  the  data  by  the  es¬ 
timate.  For  a  discussion  on  why  it  was  chosen  see  section 

4.  This  section  follows  closely  the  description  of  the  2D 
MVSE  found  in  [4]. 

The  MVSE  is  defined  as: 


Equation  10  leads  to  the  test  statistic  T(x)  = 


max.  VW  \x(m,n,t)\2 

mo, no, i  z — '  z — '  z — ' 

(m,n,t)eA(mo,no)i 


Ptgt  \ 

P mntPmnt  ) 


(ID 

The  combination  of  mo, no,  and  i  that  produces  the  largest 
test  statistic  is  the  most  likely  TGT  Model  at  the  most  likely 
initial  position  in  range  and  bearing. 

In  practice  Pmnt  is  not  known.  It  is  estimated  and  the 
data  normalized  (whitened).  After,  normalization  Pmnt  and 
Smnt  are  constant.  Therefore  , 


Ptgt  \ 
Pmnt  Pmnt  ) 


=  c  >  0 


(12) 


where  c  is  a  constant. 

So,  equation  1 1  simplifies  to, 


_Pmv(m,n,t)  — 


eHRx.'( 


(16) 


where  H  is  the  hermitian  transpose,  e  = 


.  .  z 


M  —  l  nn 


*2*3 


-444  ■ 


0  „N—1  0 


ZiZ. 


1^2 


M—l  N—l  T— In 
Z1  z2  z3  J 


with  Z\  =  exp(j2i:m/M)  ,  z^  =  exp(j2i:n/N ),  2(3  = 
exp(j2irt/T)  and  R^  is  an  estimate  of  the  autocorrelation 
matrix.  Rxx  is  composed  of  estimates  of  correlation  lags. 
In  one  dimension  the  lags  are  usually  time  so  that  the  kth 
lag  is  estimated  by  P  ^2n=o  x(n)*x(n  +  k).  In  3-D  the 
lags  are  in  lines,  beams  and  ping.  A  detailed  discussion 
of  constructing  the  R^  matrix  is  beyond  the  scope  of  this 
paper.  For  more  details  see  [4]  and  [5]. 


max 

mo, no, 


.  \x(m,n,t)\2c 

(m,n,t)eA(mo,no)i 


(13) 


Ignoring  the  constant  term  in  we  have, 


max 

mo,  no,  i 


.  \x{m,n,t)\ 2 

(m,n,t)£A(mo,no)i 


(14) 


which  is  the  test  statistic  used  in  this  paper.  Eq.  14  could  be 
written  as, 

max  Vyy/(m0,no),Wm,n,t)|2  (15) 

m0,no,i  z J  z J  z J 
( m,n,t ) 


4.  ALGORITHM 

The  algorithm  in  this  paper  is  developed  from  eq.  15.  It  is 
implemented  as  follows: 

1.  Store  N  range-bearing  plots  to  form  3-D  data. 

2.  Estimate  3-D  power  spectral  density  using  a  3-D  esti¬ 
mator  on  above  data. 

3.  Divide  data  by  the  spectral  estimate  to  generate  nor¬ 
malized  data. 

4.  Implement  the  GLRT  using  normalized  data. 


where, 

,  x  _  /  1  (m,n,t)  e  A(m0,n0)i 
I(m  o,no)i-j0  otherwise 

E  E  I(mo,  n0)i\x(m,  n,  t)\2  is  a  correlation  be- 

tween  the  data  and  the  TGT  models.  If  we  were  using  one 
TGT  model  we  would  be  maximizing  in  2-D.  The  maxi¬ 
mum  range-bearing  (mo, no)  value  would  be  the  most  likely 
initial  range  and  bearing  of  the  true  TGT.  In  this  paper  we 
are  assuming  k  TGT  models  and  therefore,  will  therefore 
maximize  over  k  sets  of  range-bearing  values. 


5.  Save  the  maximum  initial  range-bearing  cell  and  the 
associated  TGT  model.  This  is  the  range,  bearing  and 
heading  of  the  true  TGT. 

6.  Center  a  re  x  be  rectangle  over  the  maximum  range¬ 
bearing  cell  from  each  of  the  k  2-D  outputs  and  re¬ 
place  all  the  values  covered  by  the  rectangle  with  ze¬ 
ros. 

7.  Repeat  steps  5-6  until  you  have  the  desired  number  of 
possible  TGTs. 

A  more  detailed  description  is  given  next. 


4.1.  Form  3-D  Data 

To  form  3-D  data  in  a  realtime  system,  save  the  range-bearing 
plots  (pings  of  data)  until  the  desired  number  of  range-bearing 
plots  have  been  accumulated.  Then  stack  the  range-bearing 
plots  to  form  3-D  data.  In  a  non-realtime  system  (or  a  sim¬ 
ulation)  just  form  the  3-D  data  using  the  desired  number  of 
range-bearing  plots. 

4.2.  Estimate  3-D  Data 

The  data  can  be  thought  of  as  the  power  in  range,  bearing, 
and  time.  Because  of  this,  power  spectral  density  techniques 
can  be  used.  In  this  paper  the  3-D  MVSE  was  chosen,  as 
noted  above.  The  MVSE  was  chosen  because  it: 

1 .  Has  better  resolution  than  moving  average  (MA).  (see 
[4]) 

2.  Avoids  the  instability  problems  of  the  autoregressive 
(AR)  and  autoregressive  moving  average  (ARM A)  es¬ 
timators.  (see  [5]) 

3.  Is  extendable  to  any  number  of  dimensions. 

4.  Is  readily  implemented  in  MATLAB  or  C. 

An  example  of  the  spectral  estimate  is  shown  in  Figure  3, 
which  corresponds  to  the  data  in  Figure  2.  Note  Figures  2 
and  3  are  2-D  slices  of  3-D  data  and  a  3-D  spectral  estimate. 
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Fig.  2.  One  ping  slice  of  the  range-bearing-time  data  before 
normalization. 


4.3.  Generate  Normalized  Data 

Point- wise  divide  the  original  3-D  Data  by  the  estimate  of 
the  power  in  each  range-bearing-ping  cell.  The  estimate  was 
found  by  using  the  MVSE  on  the  original  data  containing 
both  target  and  noise. 
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Fig.  3.  Slice  of  the  power  spectral  density  estimate  using 
the  MVSE. 

4.4.  Correlate  Smoothed  Data 

The  GLRT  is  now  implemented  see  (15).  This  entails  corre¬ 
lating  the  normalized  data  with  the  TGT  models.  If  the  TGT 
size  (in  range  lines  and  beams),  and  the  speed  and  heading 
are  all  known,  we  would  need  only  one  TGT  Model.  In 
MATLAB  the  TGT  Model  is  a  3-D  array,  with  ones  in  the 
space  that  the  TGT  would  occupy  and  zeros  elsewhere.  For 
example,  if  the  data  set  had  four  pings,  Figure  4  is  the  TGT 
model  for  a  TGT  that  is  twenty  lines  long,  two  beam  wide, 
and  moving  up  twenty  lines  per  ping.  In  this  paper  head¬ 
ing  is  not  assumed  known.  Therefore,  we  will  need  one 
TGT  model  for  each  heading  we  wish  to  detect  a  TGT  mov¬ 
ing  in.  In  practice,  over  small  numbers  of  pings  or  for  a 
slowly  moving  TGT,  the  number  of  search  headings  would 
be  small. 

4.5.  Max  Output  of  the  Correlation 

Find  the  maximum  range-bearing  cell  over  all  the  TGT  mod¬ 
els.  This  is  the  most  likely  TGT.  This  should  have  come 
from  the  TGT  model  that  is  the  closest  match  to  the  true 
TGT  in  location  and  movement  heading.  Save  this  informa¬ 
tion  as  it  is  a  possible  TGT. 

4.6.  Zero  Out  Possible  TGT 

If  the  TGT  model  (and  the  true  TGT)  were  points  (exist¬ 
ing  in  only  one  range-bearing  cell)  then  correlating  the  data 
and  the  TGT  model  would  produce  a  largest  test  statistic 
for  the  true  TGT  with  low  outputs  elsewhere.  Since  we  are 
not  assuming  a  point  TGT  there  will  be  large  values  close 
in  range  and  bearing  to  the  true  TGT  location.  So,  if  we 
are  looking  for  the  location  of  the  N  most  likely  TGTs  we 
cannot  just  find  the  N  largest  correlation  outputs  from  each 
TGT  model,  because  one  true  TGT  or  one  area  of  high  re¬ 
verberation  will  produce  many  possible  TGTs  all  close  in 
range  and  bearing.  To  eliminate  this  possibility,  center  a 
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Fig.  4.  Example  of  TGT  model  (red  area  is  unity  and  blue 
area  is  zero) 


re  x  be  rectangle  over  the  possible  TGT  range  and  bearing 
in  each  TGT  model  and  replace  the  values  in  the  rectangle 
with  zero. 


4.7.  Repeat  Search  for  Possible  TGTs 

Eq.  15  produces  the  most  likely  true  TGT  range,  bearing 
and  heading  (TGT  model).  However,  in  practical  systems 
the  possibility  of  missing  the  true  TGT  necessitates  accept¬ 
ing  a  few  false  alarms.  Therefore,  the  N  largest  test  statis¬ 
tics  are  found  and  treated  as  possible  TGTs.  The  possible 
TGTs  could  be  passed  to  an  automatic  tracker  or  displayed 
on  a  screen. 


5.  COMPUTER  SIMULATION 


5.1.  Noise  Data 


The  data  for  each  trial  is  made  by  generating  a  3-D  array  of 
independent  non-identically  distributed  squared  CN  random 
variables.  The  size  of  the  array  is  512  range  lines,  64  beams 
and  4  pings.  The  noise  is  non-identically  distributed  be¬ 
cause  the  power  (variance)  of  the  noise  changes  with  each 
cell.  The  power  in  range  starts  low,  quickly  grows,  then 
slowly  decays  (similar  to  what  one  might  see  in  a  reverberation- 


limited  environment).  The  beams  are  combined  to  make 
range-bearing  plots  and  the  range-bearing  plots  are  com¬ 
bined  to  make  range-bearing-ping  plots  (see  Figure  1). 

5.2.  Target  Data 

The  TGT  is  then  placed  into  the  data.  The  TGT  consists 
of  the  square  of  CN  random  variables  with  power  equal  to 
twice  the  average  power  in  the  area  where  the  TGT  is  to  be 
injected.  The  TGT  is  twenty  range  lines  long  and  two  beams 
wide.  In  the  first  100  trials  the  simulated  TGT  is  moving  up 
in  beam  number  at  the  rate  of  two  beams  per  ping.  In  the 
second  100  trials  the  simulated  TGT  is  moving  up  in  line 
number  at  the  rate  of  twenty  lines  per  ping. 

5.3.  Algorithm 

The  normalization  step  is  accomplished  with  a  MVSE  as 
noted  above.  The  MVSE  used  in  the  simulation  makes  use 
of  a  one  ping,  four  lines,  and  four  beams  lag  estimate. 

The  eight  TGT  models  that  are  used  in  this  simulation 
are  matched  in  speed  to  the  injected  TGT.  They  cover  eight 
combinations  of  moving  left/right  in  beam  and  up/down  in 
line  (see  Figure  5). 


Fig.  5.  The  eight  TGT  model  vectors 

Fig  4  is  an  example  of  one  of  the  TGT  models.  This  al¬ 
gorithm  allows  for  the  addition  of  more  unknowns.  For  in¬ 
stance  if  speed  were  unknown,  a  different  TGT  set  matched 
to  each  speed  could  be  be  used.  Also  non-constant  velocity 
TGT  models  could  be  used  and  even  different  TGT  models 
for  different  size  TGTs.  But  in  this  simulation  only  the  eight 
TGT  models  mentioned  above  are  used. 

5.4.  Simulation  Details 

The  performance  of  the  algorithm  is  tested  by  repeatedly  in¬ 
jecting  a  simulated  TGT  into  the  range-bearing-ping  data  at 
a  position  when  the  data  is  reverberation  dominated.  The 
reverberation  dominated  portions  are  in  the  early  lines  (20 
to  60)  of  the  data  (see  Figure  6)  The  data  is  processed  with 
the  above  algorithm.  The  algorithm  then  produces  a  list  of 
100  possible  TGTs  sorted  by  likelihood.  The  20-line  by  2- 
beam  rectangle  (the  size  of  the  true  TGT)  is  placed  over  the 
possible  TGTs,  starting  with  the  highest  likelihood.  If  the 
true  TGT  is  within  the  rectangle  the  true  TGT  is  considered 
found.  The  possible  TGT  location  and  model  are  recorded. 
If  none  of  the  100  possible  TGTs  find  the  true  TGT  the  algo¬ 
rithm  is  considered  as  having  missed  the  true  TGT  for  that 


realization.  The  results  are  summarized  in  tabular  form  in 
Figures  7  and  8. 


caused  by  the  much  greater  changes  in  noise  power  down 
the  lines  versus  the  lesser  changes  in  noise  power  across  the 
beams. 


Lines 


6.  CONCLUSIONS 

In  this  paper  we  have  developed  an  algorithm  that  inputs 
multiple  pings  of  range-bearing  data  and  outputs  a  list  of 
possible  TGTs.  In  the  simulation  the  estimated  TGT  that 
matches  the  true  TGT  is  usually  one  of  the  highest  rank¬ 
ing  possible  TGTs.  Often  the  highest  ranking  possible  TGT 
matches  the  true  TGT.  This  means  the  true  TGT  can  be 
found  with  very  few  false  alarms.  A  complete  implemen¬ 
tation  in  MATLAB  is  available  upon  request. 


Fig.  6.  Example  of  reverberation  dominated  portion  of  one 
ping  of  data 

One  hundred  range-bearing-ping  data  sets  were  simu¬ 
lated.  The  data  sets  were  comprised  of  four  pings  worth 
of  range-bearing  plots,  each  consisting  of  512  lines  and  64 
beams. 

5.5.  Simulation  Results 

The  results  of  100  trials  where  the  TGT  is  moving  up  in 
beams  is  shown  in  Figure  7. 


true  TGT  found  in  the 

Percent 

first 

82 

2-25 

13 

26-50 

5 

51-75 

0 

76-100 

0 

Not  Found 

0 
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Fig.  7.  Performance  of  the  algorithm 

The  results  of  a  second  100  trials  where  the  TGT  is  mov¬ 
ing  up  in  lines  is  shown  in  Figure  8. 
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Fig.  8.  Performance  of  the  algorithm 

In  the  above  simulations,  the  performance  of  the  algo¬ 
rithm  for  a  TGT  moving  up  in  beams  is  better  than  the  per¬ 
formance  for  a  TGT  moving  up  in  lines.  This  could  be 


